#This code produces Table 15 in appendix F and figures 20, 21, 22, 23, 24, and 25 in Appendix G
if(!dir.exists("tabs")){dir.create("tabs")}
if(!dir.exists("tabs/appendix")){dir.create("tabs/appendix")}
if(!dir.exists("figs/appendix")){dir.create("figs/appendix")}

rm(list = ls())
library(lubridate)
library(foreign)
library(ggthemes)
library(stats)
library(lmtest)
library(sandwich)
library(quanteda)
library(stm)
library(stats)
library(ggpubr)
library(tidyverse)
library(tidytext)
gc()

# Load tishreen -----------------------------------------------------------
tishreen = read.csv("tishreen.csv")
tishreen = tishreen %>% 
  as_tibble() %>% 
  mutate(date = ymd(date),
         week = floor_date(date, unit = "week"),
         isr_pct = ofhlisrael/totalofhl) %>% 
  group_by(week) %>% 
  mutate(isr_avg = sum(ofhlisrael)/sum(totalofhl))

#Generate missing data
df = tibble(date = seq(ymd("1987-01-03"), ymd("2002-12-31"), by = "1 day"))
tishreen = left_join(df, tishreen, "date")

tishreen$year = as.factor(floor_date(tishreen$date, unit = "year"))

peace_talks = tibble(start = c(ymd("1989-12-01"), ymd("1992-07-01"), ymd("1994-04-30"), ymd("1995-05-24"), ymd("1995-12-11"), ymd("1999-12-08")),
                     end = c(ymd("1991-10-30"), ymd("1993-09-13"), ymd("1994-12-23"), ymd("1995-11-04"), ymd("1996-02-25"), ymd("2000-01-11")))

duration = c(seq.Date(peace_talks$start[1], peace_talks$end[1], by = "day"),
             seq.Date(peace_talks$start[2], peace_talks$end[2], by = "day"),
             seq.Date(peace_talks$start[3], peace_talks$end[3], by = "day"),
             seq.Date(peace_talks$start[4], peace_talks$end[4], by = "day"),
             seq.Date(peace_talks$start[5], peace_talks$end[5], by = "day"),
             seq.Date(peace_talks$start[6], peace_talks$end[6], by = "day"))
length(duration)

tishreen1 = tishreen %>% mutate(talks = ifelse(date > peace_talks$start[1] & date < peace_talks$end[1], 1,
                                               ifelse(date %in% duration, NA, 0)))
tishreen2 = tishreen %>% mutate(talks = ifelse(date > peace_talks$start[2] & date < peace_talks$end[2], 1,
                                               ifelse(date %in% duration, NA, 0)))
tishreen3 = tishreen %>% mutate(talks = ifelse(date > peace_talks$start[3] & date < peace_talks$end[3], 1,
                                               ifelse(date %in% duration, NA, 0)))
tishreen4 = tishreen %>% mutate(talks = ifelse(date > peace_talks$start[4] & date < peace_talks$end[4], 1,
                                               ifelse(date %in% duration, NA, 0)))
tishreen5 = tishreen %>% mutate(talks = ifelse(date > peace_talks$start[5] & date < peace_talks$end[5], 1,
                                               ifelse(date %in% duration, NA, 0)))
tishreen6 = tishreen %>% mutate(talks = ifelse(date > peace_talks$start[6] & date < peace_talks$end[6], 1,
                                               ifelse(date %in% duration, NA, 0)))

library(estimatr)
library(texreg)

mod1 = lm_robust(isr_pct ~ talks + year, data = tishreen1)
mod2 = lm_robust(isr_pct ~ talks + year, data = tishreen2)
mod3 = lm_robust(isr_pct ~ talks + year, data = tishreen3)
mod4 = lm_robust(isr_pct ~ talks + year, data = tishreen4)
mod5 = lm_robust(isr_pct ~ talks + year, data = tishreen5)
mod6 = lm_robust(isr_pct ~ talks + year, data = tishreen6)

texreg(list(mod1, mod2, mod3, mod4, mod5, mod6), custom.model.names = c("1989/12-1991/10", "1992/07-1993/09", "1994/04-1994/12", 
                                                                        "1995/05-1995/11", "1995/12-1996/02", "1999/12-2000/01"), 
       omit.coef = "year", custom.coef.names = c("Intercept", "Peace Talks"), float.pos = "H",
       caption = "Coverage of Israel in Tishreen during peace talks with year fixed effects", label = "tab:peace_talks_appendix", 
       fontsize = "footnotesize", include.ci = F, include.rsquared = F, include.adjrs = F, include.nobs = F, 
       reorder.coef = c(2, 1), include.rmse = F, file = "tabs/appendix/peace_talks_linear.tex")

# permutations ------------------------------------------------------------
mod1 = lm_robust(isr_pct ~ talks + year, data = tishreen1)
mod2 = lm_robust(isr_pct ~ talks + year, data = tishreen2)
mod3 = lm_robust(isr_pct ~ talks + year, data = tishreen3)
mod4 = lm_robust(isr_pct ~ talks + year, data = tishreen4)
mod5 = lm_robust(isr_pct ~ talks + year, data = tishreen5)
mod6 = lm_robust(isr_pct ~ talks + year, data = tishreen6)

beta1 = coef(lm(isr_pct ~ talks, data = tishreen1))[2]

beta2 = coef(lm(isr_pct ~ talks, data = tishreen2))[2]

beta3 = coef(lm(isr_pct ~ talks, data = tishreen3))[2]

beta4 = coef(lm(isr_pct ~ talks, data = tishreen4))[2]

beta5 = coef(lm(isr_pct ~ talks, data = tishreen5))[2]

beta6 = coef(lm(isr_pct ~ talks, data = tishreen6))[2]

nsims = 1e4
store_betas = data.frame(coef1 = numeric(nsims), coef2 = numeric(nsims), coef3 = numeric(nsims),
                         coef4 = numeric(nsims), coef5 = numeric(nsims), coef6 = numeric(nsims))

minmonth = peace_talks$start[1]
maxmonth = peace_talks$start[6]
dates = seq.Date(minmonth, maxmonth, by = "day")

s1 = length(seq.Date(peace_talks$start[1], peace_talks$end[1], by = "day")) - 1
s2 = length(seq.Date(peace_talks$start[2], peace_talks$end[2], by = "day")) - 1
s3 = length(seq.Date(peace_talks$start[3], peace_talks$end[3], by = "day")) - 1
s4 = length(seq.Date(peace_talks$start[4], peace_talks$end[4], by = "day")) - 1
s5 = length(seq.Date(peace_talks$start[5], peace_talks$end[5], by = "day")) - 1
s6 = length(seq.Date(peace_talks$start[6], peace_talks$end[6], by = "day")) - 1

set.seed(32453245)
# do the simulations
for (i in 1:nsims){
  
  # select treatment dates to kick in
  start = sample(dates, 1)
  sampled1 = seq.Date(start, start + s1, by = "day")
  sampled2 = seq.Date(start, start + s2, by = "day")
  sampled3 = seq.Date(start, start + s3, by = "day")
  sampled4 = seq.Date(start, start + s4, by = "day")
  sampled5 = seq.Date(start, start + s5, by = "day")
  sampled6 = seq.Date(start, start + s6, by = "day")
  
  d1 = tishreen %>% mutate(treated = ifelse(date %in% sampled1, 1, ifelse(date %in% duration, NA, 0)))
  d2 = tishreen %>% mutate(treated = ifelse(date %in% sampled2, 1, ifelse(date %in% duration, NA, 0)))
  d3 = tishreen %>% mutate(treated = ifelse(date %in% sampled3, 1, ifelse(date %in% duration, NA, 0)))
  d4 = tishreen %>% mutate(treated = ifelse(date %in% sampled4, 1, ifelse(date %in% duration, NA, 0)))
  d5 = tishreen %>% mutate(treated = ifelse(date %in% sampled5, 1, ifelse(date %in% duration, NA, 0)))
  d6 = tishreen %>% mutate(treated = ifelse(date %in% sampled6, 1, ifelse(date %in% duration, NA, 0)))
  
  # run regression
  mod1 = lm(isr_pct ~ treated, data = d1)
  mod2 = lm(isr_pct ~ treated, data = d2)
  mod3 = lm(isr_pct ~ treated, data = d3)
  mod4 = lm(isr_pct ~ treated, data = d4)
  mod5 = lm(isr_pct ~ treated, data = d5)
  mod6 = lm(isr_pct ~ treated, data = d6)
  
  store_betas$coef1[i] = coef(mod1)["treated"]
  store_betas$coef2[i] = coef(mod2)["treated"]
  store_betas$coef3[i] = coef(mod3)["treated"]
  store_betas$coef4[i] = coef(mod4)["treated"]
  store_betas$coef5[i] = coef(mod5)["treated"]
  store_betas$coef6[i] = coef(mod6)["treated"]
}


# get one-sided p-value 
pval1 = mean(store_betas$coef1 < beta1, na.rm = T)
pval2 = mean(store_betas$coef2 < beta2, na.rm = T)
pval3 = mean(store_betas$coef3 < beta3, na.rm = T)
pval4 = mean(store_betas$coef4 < beta4, na.rm = T)
pval5 = mean(store_betas$coef5 < beta5, na.rm = T)
pval6 = mean(store_betas$coef6 < beta6, na.rm = T)

# tishreen plot results ------------------------------------------------------------
p1 = ggplot(store_betas) + 
  aes(x = coef1) + 
  geom_histogram(fill = "#A0A0A0", aes(y = ..ncount../sum(..ncount..)), bins = 40) + 
  geom_segment(x = beta1, xend = beta1, y = 0, yend = .07, size = 1) + 
  annotate(geom = "text", x = beta1, y = .08, label = paste0("1989-1991\nOne-sided\np = ", round(pval1, 3))) + 
  labs(x = "Placebo coefficient estimate", y = "Density") +
  scale_x_continuous(breaks = seq(-.1, .1, .05)) + 
  theme_minimal() + 
  theme(panel.grid = element_blank(),
        panel.border = element_rect(colour = "black", fill=NA, size=.85),
        panel.spacing = unit(0, "in"))
p1
  ggsave("figs/appendix/tishreen_placebo1.pdf", width=11, height = 7)

p2 = ggplot(store_betas) + 
  aes(x = coef2) + 
  geom_histogram(fill = "#A0A0A0", aes(y = ..ncount../sum(..ncount..)), bins = 40) + 
  geom_segment(x = beta2, xend = beta2, y = 0, yend = .07, size = 1) + 
  annotate(geom = "text", x = beta2 - 0.02, y = .074, label = paste0("1992-1993\nOne-sided\np = ", round(pval2, 3))) + 
  labs(x = "Placebo coefficient estimate", y = "Density") +
  scale_x_continuous(breaks = seq(-.1, .1, .05)) + 
  theme_minimal() + 
  theme(panel.grid = element_blank(),
        panel.border = element_rect(colour = "black", fill=NA, size=.85),
        panel.spacing = unit(0, "in"))
p2
ggsave("figs/appendix/tishreen_placebo2.pdf", width=11, height = 7)


p3 = ggplot(store_betas) + 
  aes(x = coef3) + 
  geom_histogram(fill = "#A0A0A0", aes(y = ..ncount../sum(..ncount..)), bins = 40) + 
  geom_segment(x = beta3, xend = beta3, y = 0, yend = .07, size = 1) + 
  annotate(geom = "text", x = beta3 - 0.02, y = .074, label = paste0("1994-1994\nOne-sided\np = ", round(pval3, 3))) + 
  labs(x = "Placebo coefficient estimate", y = "Density") +
  scale_x_continuous(breaks = seq(-.1, .1, .05)) + 
  theme_minimal() + 
  theme(panel.grid = element_blank(),
        panel.border = element_rect(colour = "black", fill=NA, size=.85),
        panel.spacing = unit(0, "in"))
p3
ggsave("figs/appendix/tishreen_placebo3.pdf", width=11, height = 7)

p4 = ggplot(store_betas) + 
  aes(x = coef4) + 
  geom_histogram(fill = "#A0A0A0", aes(y = ..ncount../sum(..ncount..)), bins = 40) + 
  geom_segment(x = beta4, xend = beta4, y = 0, yend = .07, size = 1) + 
  annotate(geom = "text", x = beta4 - 0.02, y = .074, label = paste0("1995-1995\nOne-sided\np = ", round(pval4, 3))) + 
  labs(x = "Placebo coefficient estimate", y = "Density") +
  scale_x_continuous(breaks = seq(-.1, .1, .05)) + 
  theme_minimal() + 
  theme(panel.grid = element_blank(),
        panel.border = element_rect(colour = "black", fill=NA, size=.85),
        panel.spacing = unit(0, "in"))
p4
ggsave("figs/appendix/tishreen_placebo4.pdf", width=11, height = 7)

p5 = ggplot(store_betas) + 
  aes(x = coef5) + 
  geom_histogram(fill = "#A0A0A0", aes(y = ..ncount../sum(..ncount..)), bins = 40) + 
  geom_segment(x = beta5, xend = beta5, y = 0, yend = .07, size = 1) + 
  annotate(geom = "text", x = beta5 - 0.03, y = .074, label = paste0("1995-1996\nOne-sided\np = ", round(pval5, 3))) + 
  labs(x = "Placebo coefficient estimate", y = "Density") +
  scale_x_continuous(breaks = seq(-.1, .1, .05)) + 
  theme_minimal() + 
  theme(panel.grid = element_blank(),
        panel.border = element_rect(colour = "black", fill=NA, size=.85),
        panel.spacing = unit(0, "in"))
p5
ggsave("figs/appendix/tishreen_placebo5.pdf", width=11, height = 7)

p6 = ggplot(store_betas) + 
  aes(x = coef6) + 
  geom_histogram(fill = "#A0A0A0", aes(y = ..ncount../sum(..ncount..)), bins = 40) + 
  geom_segment(x = beta6, xend = beta6, y = 0, yend = .07, size = 1) + 
  annotate(geom = "text", x = beta6 - 0.03, y = .074, label = paste0("1999-2000\nOne-sided\np = ", round(pval6, 3))) + 
  labs(x = "Placebo coefficient estimate", y = "Density") +
  scale_x_continuous(breaks = seq(-.1, .1, .05)) + 
  theme_minimal() + 
  theme(panel.grid = element_blank(),
        panel.border = element_rect(colour = "black", fill=NA, size=.85),
        panel.spacing = unit(0, "in"))
p6
ggsave("figs/appendix/tishreen_placebo6.pdf", width=11, height = 7)



